Nonlinear evolution of the step meandering instability of a growing crystal surface 
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The growth of crystal surfaces, under non-equilibrium conditions, involves the displacement of 
mono-atomic steps by atom diffusion and atom incorporations into steps. The time-evolution 
of the growing crystal surface is thus governed by a free boundary value problem [known as the 
Burton-Cabrera-Franck model]. In the presence of an asymmetry of the kinetic coefficients [Erlich- 
Schwoebel barriers], ruling the rates of incorporation of atoms at each step, it has been shown that 
a train of straight steps is unstable to two dimensional transverse perturbations. This instability is 
now known as the Bales-Zangwill instability (meandering instability). We study the non-linear evo- 
lution of the step meandering instability that occurs on a crystalline vicinal surface under growth, 
in the absence of evaporation, in the limit of a weak asymmetry of atom incorporation at the steps. 
We derive a nonlinear amplitude equation displaying spatiotemporal coarsening. We characterize 
the self-similar solutions of this equation. 

PACS numbers: 81.15.Hi, 68.35.Ct, 81.10.Aj, 47.20.Hw 
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Molecular beam epitaxy (MBE) is often used to grow 
nano-structures on vicinal surfaces of semiconductor and 
metallic crystalsi2i2iii^&. The ability to grow smooth 
crystals with a sharp interface is of considerable impor- 
tance when manufacturing electronic and opto-electronic 
devices. One possibility for achieving this goal is to use 
the step-flow mode, where deposited adatoms diffuse and 
attach directly to preexisting steps on a vicinal surface. 
Ideally the surface grows without changing its shape by 
advancement of a uniform train of step traveling at con- 
stant velocity. In experiments, a wide number of addi- 
tional effects may modify this scenario. Fluctuations in 
the beam intensity and island growth may lead to kinetic 
roughening of the surface^. The presence of impurities 
can also pin the steps at random positions. More impor- 
tantly, intrinsic morphological instabilities such as step 
bunching and step meandering can take place on a vici- 
nal surface and lead to the destruction of a stable train 
of equidistant straight steps. Thus, under standard MBE 
growth conditions (one monolayer per minute) a rich va- 
riety of crystal surface morphologies is experimentally 
observed, resulting from the nonlinear evolution of step 
bunching and meandering instabilities£*2i2iiS. Moreover, 
the self-organization arising from these instabilities has 
been proposed as a natural candidate for the develop- 
ment of technological applications such as quantum dots 
and quantum wells 1112 . 

The step meandering instability was originally pre- 
dicted theoretically by Bales and Zangwill 13 for a vic- 
inal surface under growth. Its origin comes from the 
asymmetry between the descending and ascending cur- 
rents of adatoms. As shown by Bales and Zangwill, a 
straight train of steps during MBE growth may become 
morphologically unstable in the presence of a kinetic at- 
tachment asymmetry at the step: the Erlich-Schwoebel 
effect (ES). The physical mechanisms taking place during 
this instability are the destabilizing ES effect combined 
with the stabilizing effect provided by the step rigidity. 



It was shown that the most dangerous mode is the syn- 
chronized mode for which all the steps have the same 
phased. Nonlinear extensions of this work have shown 
that the meander evolution can be described by ampli- 
tude equations showing various behaviors. Close to the 
instability threshold, starting from the Burton-Cabrera- 
Frank (BCF) model, it was showni^ that the step po- 
sition in the presence of desorption (evaporation) obeys 
the Kuramoto-Sivashinsky equation: 



(dyhf 



(1) 



where y is the coordinate along the step and x = h(y, t) 
is the meander amplitude which describes the step shape 
in the (x,y) plane. The ultimate stage of this dynam- 
ics is thus found to be spatiotemporal chaos. In the 
case of negligible desorption with strong or moderate ES 
effecU£ilU&, it was shown that the step shape obeys a 
highly nonlinear equation 



d t h = -d,. 
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This equation cannot be derived from a weakly nonlinear 
analysis but is based on the assumption that the slope 
of the steps is order unity. Instead of spatiotemporal 
chaos, a regular pattern is revealed, the lateral mod- 
ulation wavelength is fixed while the amplitude of the 
step deformation (transverse meandering amplitude) in- 
creases like t 1 / 2 . Later, it was shown that the inclusion 
of the elastic step interactions affects the step dynam- 
ics in the sense that they induced a lateral coarsening 19 . 
Finally, it was recently shown that interrupted coarsen- 
ing occurs when two-dimensional diffusion anisotropy is 
included^!. 

In this paper, we show using standard weakly nonlin- 
ear analysis, that under the assumption of negligible des- 
orption and weak Erlich-Schwoebel effect (ES), the time- 
evolution of the Bales-Zangwill instability is governed by 
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the following equation: 



d t h = —d y h — d y h 



d 2 y {d y hf 



(3) 



This equation was already mentioned in Ref. 18 on the ba- 
sis of symmetry arguments as a possible candidates for 
the time evolution of the meandering amplitude but was 
not explicitly obtained due to a different choice of the 
physical parameters (large and order one ES effects)^. 
Our results are illustrated by numerical simulations. We 
show using a simple similarity argument, that the char- 
acteristic transverse meandering amplitude (step width) 
grows like t and that the characteristic lateral coarsening 
exponent grows like i 1 / 2 . Our result are in agreement 
with the numerical simulations performed in Ref^. 

Let us denote by x n (y, t) the positions at time t of the 
n step (see Fig. 1). For simplicity we neglect elastic in- 
teractions between steps and assume that the desorption 
of adatoms and transparency of the steps are negligi- 
ble. During growth, the adatom surface concentration 
on each terrace C n (x, y, t) , obeys the following diffusion 
equation!^: 



DV 2 C n (x,y)+F = 0, 



(4) 



where D is the adatom diffusion coefficient, and F the 
deposition flux. This equation for C n is supplemented by 
the following boundary conditions: 

Dfl ■ VC„ = V + {C n - Ce«,n-l) , X = X n -l{V,t) , (5) 
Dn-WC n = -V-(C n - C eq ,n) , X = X n (y,t) . (6) 

Here v + and v _ are the ES coefficients which are propor- 
tional to the rate of attachment of adatoms on the steps 
from the terrace; fi is the external normal to the step : 



(1, -dyX n )/J 1 + (dyXn) 2 , 



(7) 



and V is the two dimensional gradient operator (d x , d y ). 
The adatom equilibrium concentrations C eq . n depend on 
the step curvatures KrP^'- 



Ceq.n = C(0) (1 + ^ K n) , 



(8) 



where T = flj/ksT, with Q the unit atomic surface, 
T the temperature, ks the Boltzmann constant, 7 the 
step rigidity and C(o) the reference adatom concentra- 
tion. The step curvature is given by 



d 2 x 



(9) 



(1 + {d y x n yf/ 2 • 

The normal velocity of the n step is given by 

v n = DVL(ii- VC n+ i -n- VC„) , x = x n (y,t). (10) 
The x component of the normal velocities reads: 

X n =V n (l+ (dyXr,) 2 ) 1 ' 2 . (11) 




FIG. 1: Sketch of the vicinal surface showing a succession of 
steps and terraces and the development of the meandering 
instability. x n and ^ n (j/,i) are step positions and the cor- 
responding in-phase perturbations. C„ denotes the adatom 
concentration on the terrace. 



where the dot represents the time derivative. 

In order to get a non-dimensional version of equations 
(|4I5I8I10|I . we set the unit of length to be the initial size 
of the terrace Iq (initial distance between steps) and the 
unit of time to be Iq/ (C(o)FttD). Specifically, one sets 



~ - x ~ - y r - MC» ~ C(o)] 

Iq Iq ^O) 1 



t = 



mc {0) Dr 



(12) 

in equations if^-lf^-ll^- lfTrHl . We obtain after omitting 
the tildes on the variables the following equations for the 
dimcnsionless variables 

V 2 C„ = -/ (13) 

n-VC„ = a + (C n - k„-i) ,x = x n -i(y,t) (14) 

n-VC„ = -a_(C„-K„) , x = x n (y,t) (15) 

v n = h ■ VC„+i - h ■ VC„ , x = x n (y, t) (16) 

The system is thus controlled by three independent pos- 
itive nondimensional parameters: 



f 



Fl 3 



C {0) TD 



D 



a.- 



D 



(17) 



respectively related to the flux and attachment lengths. 
Let us investigate the linear stability of a train of equidis- 
tant steps traveling at a constant velocity / when per- 
turbed transversally. The shape of the steps can be de- 
composed in Fourier modes of the form 



x n (y,t) =n + ft + £ n (y,t) 



where 



Cn(*,y) 



— „<?(q,<f>)t+iqy+in</> 



(18) 



(19) 



where q and <fi & r e respectively the wavenumber and the 
phase of the perturbation (see Fig. 1). Inserting these 
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expressions into equations (|13ll4ll5ll6fl . we obtain the 
general dispersion relation a = a(q, <f>) which possesses, 
for each (j>, one maximumiii^. The maximum growth 
rate is reached for the in-phase perturbation <j> — 0. Al- 
though the full expression of the dispersion relation is 
cumbersome, near the instability threshold ( / = 0), we 
introduce a small parameter e measuring the distance to 
the threshold and we also assume that the E-S effect is 
small. This latter assumptions was not used in previous 
worksi&I&. The small parameter e arises naturally when 
considering the long wavelength limit, in which q — ► eq. 
We therefore set 

f = ef, C = e 2 C. a + = a_+e 3 5, t = 4 (20) 

in equations I|13I14I14I16[) . This scaling will lead to the 
following equations after omitting the tilde 

V 2 C„ - -fe (21) 

n-VC„ = a + (C n - e 2 K„_i) , x = x n -i(y, t) (22) 

n ■ VC„ = -a-(C n - e 2 K n ) , x = x„(y,t) (23) 

v n = n ■ VC„ + i - n ■ VC„ , x = x n (y, t) .(24) 

To lowest order in e we find that the growth rate a = 
a(q,0) is, 

a = e 6 (W(a_+2)^-/) , (25) 

The growth rate a is maximum for the wavenumber 
<Zmax = V fSa-(2 + a_)/2, so that the most unstable 



d t u = -d 2 



where a = a_. We have renamed the capital letters T 
and 7 as t and y for simplicity. It is necessary to go to 
the seventh order to balance the nonlinear term with the 
linear ones. This equation was recently presented in the 
context of surface growth on a vicinal surface of silicon 
characterized by different types of steps and terraces^ 5 -. 
In this previous work, it was shown that the anisotropy 
of diffusion induces a meandering instability and sim- 
ple similarity and matching arguments lead to a com- 
plete picture for the long time behavior of the solution of 
equation (|27|l . It is worth noting that Eq. I|27|) admits 
an exact particular solution in the form of a stationary 
parabola 

, , 3a 2 (2 + a) 2 S , 

u »' f = ~ Affile 2 2^ ■ 28 
4(6 -I- ba + cr) 

Trying a similarity solution u(y,t) = t a ip(y/t b ) of equa- 



wavelength A should scales ~ \f\f]. The meandering 
instability originates from the asymmetry between the 
descending and ascending currents of adatoms; the in- 
stability occurs when the rate of attachment of adatoms 
from the terrace lower a+ is greater than the rate of 
attachment of adatoms from the upper terrace. In the 
opposite case, the meandering instability is not present 
but a step bunching instability develops. However, at 
this moment due to the technical difficulties of the ex- 
periments there has not yet been a clear demonstration 
of the Bales-Zangwill instability. 

We study now the nonlinear evolution of the mean- 
dering instability, in the limit of weak amplitudes and 
long wavelengths. In order to obtain the relevant non- 
linear dynamics we use a standard multi-scale method. 
Adatoms concentrations and step shapes are expanded 
in powers of e. The adatoms concentrations depends on 
the slowly varying space and time variables: Y = ey and 
T = e 4 t. The details of the calculations are deferred to 
Appendix A. We thus write 

7 

C n {x,y) = Y, e%C n\ x ' Y ) ,Uv,t) = eu(Y,T). (26) 

i=l 

Solving diffusion equation (|21() and boundary conditions 
(|22I23|) up to order e 7 , and inserting the results into the 
step velocity equations l|24|l , we find the equation for the 
slowly varying amplitude at order e 7 . It may be written 
as : 



(27) 

I 

tion (|27fl . we obtained the exponents = 1/2 and a = 1, 
which agree with our numerical results 2 ^. The general, 
asymptotic solution of equation l|27|l can be thought of 
a superposition of parabolas as shown on the Figure 3 
of Ref»2^. The joining regions between parabolas can be 
matched using a solution of the Burgers equation ob- 
tained by a Hopf-Cole transformation 25 . A typical spa- 
tiotemporal evolution from a random initial condition is 
shown in Fig. 2. The step width (transverse meandering 
amplitude) of the steps increase linearly in timeS 5 . and 
this results is in agreement with the results obtained by 
Kinetic Monte-Carlo methods in the Fig. 7 of RefiSS. 

In this article, we have shown that the amplitude of 
the meanders is governed by Eq. (|27|l which displays 
non-interrupted coarsening. It would be interesting in 
a further study to investigate the effect of phase free- 
dom of the steps and to derive a continuous coarse- 



'Sf 



a(a + 2)u + d 2 u + 



f[a(a + 6) +6] 
6a(a + 2) 



{d y uf 
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FIG. 2: Spacetime plot of u(y, t) with nondimensional y and 
t axes, given by the numerical solution of equation 1271 with 
a = 1, S = 1 and f — 1. The coarse-graining of structures 
leads to a superposition of parabolas"*, with a size (u 2 } 1//2 ~ 
t. In the long-time state all the parabolas tend to have unity 
curvature at their maximum, and width increasing as \ft. 

grained model for the space-time evolution of a vicinal 
surface. Such coarse-grained models could also be ap- 
plied to step bunching 2 ^£ and to electromigration in- 
duced instabilities^^. We hope that the present work 
will motivate more experimental research in the subject 
and we plan to extent this work to the study of the 
coupling between step bunching instabilities and step 
meandering instabilities. It would also be interesting 
to study, using a similar line of thought, the effect of 
an elastic stress commonly encountered in heteroepitaxy 
experiments^*^ under step flow conditions^. It would 
also be interesting to find out if there is a cross-over be- 
tween the weak nonlinear regime we have described in 



this paper and the strongly nonlinear regime described 
in RefsiSiiLia. 
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APPENDIX A: DETAILS OF THE EXPANSION 

From equations (|21I22I23I26[) . we obtain the following 
relations: 

d 2 x CU = -f, (AI) 

d 2 x C&=0, (A2) 
C*(x,Y) = ±fu(Y)(2x + l), 



%CW+%CP=0, (A3) 



^+W=0, (A4) 



Ck 4) (x,Y) = 



6/Qa - 1)6 - f(2x + l)a 2 [x(x + l){a + 2) - 1] u"(Y) 
I2a 2 (a + 2) 
I 



(A5) 



cL 5 Hx,Y) = 



3(a + 2) [f(2x(x + l)a - l)u'(Y) 2 - 4au"(Y)] + fu [a [a + 6x(x + l)(a + 2)]u"(Y) 



12a(a 
I 



2) 



(A6) 
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C<?\x,Y) = 



f(2x + 1) (360u(Yy 2 (Y){a + 2) 2 + 180u 2 {Y)u"(Y)(a + 2) 2 ) 

720(a + 2) 2 

/(2s + 1) [- (a + afo + l)(g + 2)(-a + 3x(ir + l)(g + 2) - 12) + 12) u""(y)]) 

720(a + 2) 2 



r 



(A7) 



cl r) (*,y) = 



(720a 3 (a + 2) 2 ) 



(A8) 



7 = 360/(a + 2) 2 u 2 (r)u' 2 (F)a 3 + 120/(a + 2) 2 u 3 (y)u"(r)a 3 

-/ (30(a + 2) 2 .t 4 + 60(a + 2) 2 .t 3 + 30a(a + 2)a; 2 - 60(a + 2)x - a(a + 12)) w(y)u""(y)a 3 
-30 (/a(a + 2)(a(2(a + 6) + x(x + l)(a(3a;(a; + l)(a + 2) - 2 (a + 6)) - 12)) + 12)m" 2 (T)) 
-30 (-2/a (a (3(a + 2)a; 2 + 2(a + 3)i - 2) - 6) fe"(F)) 
-30 (2 (6/(a + l)(xa- 1)<5 2 )) 

-30 (+a(a + 2) (/(a(a + i(a; + l)(a(-a + 2x(a; + l)(a + 2) - 7) - 6) + 6) + 6)u'(y)u"'(Y))) 
-30 (-6a(a + 2){x{x + l)a - l)«""(y)) , 

I 



Here the symbol ' denotes the derivative with respect 
to the variable Y. 
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